Characterization of the Melting Transition in Two Dimensions 
at Vanishing External Pressure Using Molecular Dynamics Simulations 
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A molecular dynamics study of a two dimensional system of particles interacting through a 
Lennard- Jones pairwise potential is performed at fixed temperature and vanishing external pressure. 
As the temperature is increased, a solid-to-liquid transition occurs. When the melting temperature 
T c is approached from below, there is a proliferation of dislocation pairs and the elastic constant 
approaches the value predicted by the KTHNY theory. In addition, as T c is approached from above, 
the relaxation time increases, consistent with an approach to criticality. However, simulations fail 
to produce a stable hexatic phase using systems with up to 90,000 particles. A significant jump in 
enthalpy at T c is observed, consistent with either a first order or a continuous transition. The role 
of external pressure is discussed. 



I. INTRODUCTION 

Melting of an infinite solid in two dimensions has 
been described as a process driven by a proliferation 
of thermally excited dislocation pairs in the Kosterlitz- 
Thouless-Halperin-Nelson and Young (KTHNY) the- 
orySHi. The theory, formulated at vanishing external 
pressure, predicts the existence of a new, "hexatic" , inter- 
mediate thermodynamic phase. While solids are charac- 
terized by long range translational and orientational or- 
der, liquids only present short range order. The predicted 
hexatic phase presents long range orientational order but 
lacks long range translational order. The KTHNY theory 
predicts a second order phase transition from the crys- 
talline to the hexatic phase at which point there is a uni- 
versal jump of a normalized elastic constant from a finite 
value to zero. This first transition is followed by a second 
transition from the hexatic phase to the liquid phase at 
a higher temperature. The KTHNY theory continues to 
generate interest, specially because increased numerical 
capabilities and new experimental techniques currently 
allow for new and more accurate testing of theoretical 
predictions. Indeed, there have been numerous attempts 
at the verification, both experimentally and numerically, 
of the KTHNY theoretical predictions, with mixed out- 
comes. 

On the experimental side, studies with colloidal par- 
ticles have provided evidence of two stage melting with 
an intermediate hexatic phase^^l, and of elasticity be- 
havior in agreement with the KTHNY predictions®. 
The obse rved transition, however, appears to be first 
ordeJSlLQ]. Similar results have been obtained with diblock 
copolymers^. Recently, melting in two steps with an in- 
termediate hexatic phase has been observed in monolay- 
ers of polycristalline colloidal films, but not in thin or 
thick multilayer films^. Also recently, but in a differ- 
ent context, dislocations have been directly observed in 
graphene^, prompting a renewed interest on the role of 
defects in this two dimensional materiaP^^. 



On the numeric al sid e, molecular dynamics and Monte- 
Carlo simulation ;] 19 * 20 ! of systems with a small number 
of particles (N) broadly detected a transition where the 
number of dislocations proliferates, but failed to pro- 
vide clear evidence for the nature of the observed tran- 
sition. F irst order melting has been reported in the 
literature) 19 * 21 * ^, whil e other calculations support a con- 
tinuous transitio n 1 24 * 25 ^. 

The critical properties of the KTHNY transition are 
a consequence of the renormalization effect that small 
scale fluctuations have on large scale fluctuations. For 
this mechanism to be operative, well separated length 
scales must exist, suggesting a minimum size for numer- 
ical simulations in two dimensions of 10 . Indeed, Chen 
et al. ^ performed molecular dynamics simulations of a 
Lennard- Jones system with a varying number N of par- 
ticles. They found a metastable hexatic phase for sys- 
tems with N > 36, 864, but not for N < 16, 386. In all 
cases, simulations were performed at a significant exter- 
nal pressure, a fact that alters the dislocation generation 
mechanism: the interaction between the components of a 
dislocation pair tends to close it down, while the external 
pressure, for some orientations, tends to open it up. The 
whole process becomes one of thermal activation, much 
like nucleation, and the likelihood of having isolated 
dislocations — and an hexatic phase — increases. A subse- 
quent study in terms of inherent structure theory showed 
consistency with the KTHNY theory 2 ^. More recently, 
a molecular dynamics study 2 ^* carried out at constant 
volume, involving 36,000 particles interacting through a 
Lennard- Jones potential also reported the presence of an 
hexatic phase between the solid and liquid phases. How- 
ever, phase coexistence in NVT ensembles precludes un- 
ambiguous interpretation of these results. 

On a different vein, in a three dimensional continuum 
elastic solid, dislocation lo ops d rive a mechanical instabil- 
ity at a finite temperatur e 29 * 30 !, at which point the shear 
modulus vanishes as a function of reduced temperature, 
following a power law with an exponent whose value is a 
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function that is independent of the microscopic details 
of the elastic solid. Numerical calculations of super- 
heated Lcnnard- Jones crystals near the melting transi- 
tion in three dimensions show the a ppear ance of dislo- 
cations as the temperature is raisecpIEZl. However, it 
is unclear whether these dislocations play a central role 
during the phase transition or are just a by-product of 
another transition-driving mechanism. 

The present work used molecular dynamics (MD) sim- 
ulations in the NpT ensemble to study the melting tran- 
sition in two dimensions. MD has been chosen in this 
occasion over the possible Monte Carlo alternative in 
order to study the time evolution of the system, espe- 
cially its relaxation behavior near criticality. We found 
a single-step solid-to-liquid transition (as determined by 
the enthalpy change) when a vanishing external pressure 
was applied to the system, in contrast to a multi-step 
transition at high pressure like the one presented irP^ 
using the same number of particles. However, within a 
narrow temperature interval defining the solid-to-liquid 
transition, the monitored relaxation times, elastic con- 
stants, and the evolution in number of dislocations were 
all consistent with the KTHNY theory. We suggest that, 
because of the necessary interplay between many length 
scales, a stable hexatic phase will be unambiguously ob- 
served only in systems with at least ~ 10 6 particles at 
zero external pressure. 



0.005% per million iterations. We also verified momen- 
tum conservation and monitored pressure and tempera- 
ture throughout simulations. Thermostat and barostat 
frequencies oj p and uj$^ were set to values between 80 
and 100, and 0.1 and 0.4, respectively. The cutoff radius 
was set to 4 and a Verlet-neighbours list was used with 
a radius of 5.9. 

Initial conditions were set by initial positions corre- 
sponding to a perfect triangular lattice, or a perfect tri- 
angular lattice plus a randomly oriented displacement 
of 0.05 LJ units, as well as initial velocities given by a 
Gaussian distribution corresponding to each temperature 
(initial condition IC1). Otherwise, initial conditions were 
assigned by equilibrium values of positions and velocities 
obtained from a previous run with similar values of T 
and p (initial condition IC2). Most runs were carried out 
with TV = 36,864, with up to 1.7 x 10 7 time iterations, 
and a smaller number of runs with N = 90, 000. 

To characterize the state of the system throughout sim- 
ulations we monitored the time evolution of the system's 
enthalpy and computed the pair and orientational cor- 
relation functions g(r) and ge(r), respectively^. The 
Lame parameters were also computed, as a function of 
temperature. The number of dislocation pairs present in 
the crystal below the transition was examined counting 
the number of nearest neighbors for each particle, and 
through visualization with the aid of a Voronoi construc- 
tion. 



II. MOLECULAR DYNAMICS SIMULATIONS 

Molecular dynamics (MD) simulations were carried out 
using a parallel MD code developed "in house" based on 
the libraries presented irP^. Simulation systems com- 
prised N identical particles of mass m in two dimensions 
interacting through a pairwise truncated and shifted 
Lennard- Jones (LJ) potential, 



</>ljM = 4e 



12 



-c 



(1) 



for inter-particle distances r smaller than a cutoff ra- 
dius r c and for r > r c . The value of C was chosen 
to ensure continuity of the potential. All simulations 
were performed using periodic boundary conditions with 
a fully-flexible cell in the NpT ensemble defined by the 
non-Hamiltonian equations of motion described by Mar- 
tyna, Tobias, and KleirPS. The equations of motion were 
integrated using the 5-value Gear predictor-corrector al- 
gorithm. 

All simulation parameters and monitored quantities 
are expressed in reduced units: x = x*a; t = t* \J ma 2 Je; 
T = T*e/kB, where fcs is the Boltzmann constant; E — 
E*e; and p = p*e/a 3 (in the case of Argon, e = 0.0104eV 
and a = 3.4 A, leading to a a unit of temperature T be- 
ing 120.6 K and a unit of pressure p being 42 MPa). The 
integration time step was set to O.OOOS-^/mo^/e (34 fs 
for Argon)^, ensuring extended-energy conservation to 



III. RESULTS 

A. Simulations at different temperatures show 
melting transition 

Two sets of simulations were performed to study the 
melting behavior of a LJ system. The first set of simu- 
lations, aimed at reproducing the results presented irP^, 
were performed at p — 20 with N = 36, 864. Three sim- 
ulations performed at temperatures T\ — 2.15, T^ = 2.16 
and T 3 = 2.17 were carried out using a perfect crys- 
talline lattice as initial condition. Despite the use of a 
different set of equations of motion that are modularly 
invariant, we obtained results that are consistent with 
those presented irP^. The enthalpy h as a function of 
time remained stable during the simulation at T\, in- 
creased rapidly in an apparent single step to achieve a 
stable value during the simulation performed at T 3 , and 
increased in two steps with a transient state in the sim- 
ulation performed at temperature Ti (data not shown). 
The pair and orientational correlations functions, g(r) 
and gs(r), were consistent with a solid phase for the 
simulation performed at T\ and with a liquid phase at 
the end of simulations performed at T2 and T3. The 
transient state observed at T2 exhibits long-range orien- 
tational but not translational order, consistent with an 
hexatic phase. Similar results were obtained with a simu- 
lation performed at T2 = 2.16 that used a thermalized (at 
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FIG. 1: Enthalpy as a function of time for several temper- 
atures. Upper panel: initial conditions are given by initial 
positions corresponding to a perfect triangular lattice plus a 
randomly oriented displacement of 0,05 in LJ units, and initial 
velocities given by a Gaussian distribution corresponding to 
each temperature (initial condition IC1). From the bottom- 
right to the top-left curves the temperatures are T=0.4050, 
0.4095, 0.4130, 0.4160, 0.4200, and 0.4250. Lower panel: posi- 
tions and velocities are provided, at each temperature, by the 
equilibrium values obtained by a previous run at T — 0.40725 
(initial condition IC2). From the bottom-right to the top- 
left curves the temperatures are T=0.40725, 0.41000, 0.41250, 
0.41500, 0.42000, 0.42500, and 0.44000. All simulations were 
carried out at vanishing external pressure. 



T = 2.15) set of initial conditions. These results show the 
possible existence of a metastable and transient hexatic 
phase when simulations are performed at a high external 
pressure p = 20, as presented irPH 

To test the behavior of the system at a vanishing exter- 
nal pressure, we performed a second set of simulations in 
which a system with 7Y = 36, 864 particles starting from 
initial conditions IC1 and IC2 was simulated at several 
different temperatures and pressure p = (more pre- 
cisely, with vanishing normal and tangential stresses). 
As with the first set of simulations, the enthalpy h of 
the system as a function of time (Figure 1) was sta- 



ble at low temperature (T < 0.40725), but increased 
in an apparent single-step transition to an equilibrium 
value for high temperatures (T > 0.4095). The nature 
of the phase was characterized using g(r) and g§{r) at 
the ending configurations of each simulation, confirm- 
ing a liquid phase for T > 0.4095, and a solid phase 
for T < 0.40725 (Figure [2]). As opposed to the results 
obtained at p — 20, no intermediate and transient state 
(possibly corresponding to a hexatic phase) was observed. 
It is possible that a hexatic phase with an algebraically 
decaying ge(r) and an exponentially decaying g(r) could 
be found within the interval of temperatures given by 
T_ = 0.40725 < T < T + = 0.4095. However, it was 
not possible to reach equilibrium in between these two 
temperature values within the available simulation time 
scales (data not shown). Defining the relative change 
in enthalpy as Ah = 2(h + — h-)/(h + + (where 
h± = h(T±)) and the relative change in temperature as 
AT = 2(T + - T-)/(T+ + T_), we find Ah = 0.1433 and 
AT = 0.0044. Within this narrow temperature interval, 
our data is consistent both with an abrupt jump from a 
low to a high enthalpy value at some intermediate tem- 
perature, as would be the case for a first-order transition 
(i.e., with latent heat), as well as with a two-step change 
of enthalpy as a function of temperature within the tem- 
perature interval, as would be the case for a continuous 
transition, without latent heat (Figure (31). 



B. Relaxation times slightly above melting increase 
as the melting temperature is approached 

Near a critical point, relaxation times increase as crit- 
icality is approached. This is due to the increasing size 
of fluctuations^], that reach macroscopic dimensions at 
the critical point. To determine the nature of the phase 
transition observed at p = 0, we monitored the relaxation 
time (tpi) as a function of temperature in the simulations 
mentioned at T > T c (Figure [2]). 

The relaxation time tn is defined, grossly, as the in- 
stant where the second derivative of enthalpy vs time 
vanishes, and more precisely as follows: the curve en- 
thalpy vs time is interpolated by a smooth function whose 
second derivative is computed numerically, and the times 
i mEtx , where curvature is a maximum and t min , where it 
is a minimum, determined. The relaxation time is then 
defined through t R = i max + (t min - t max )/2. 

Within the accuracy of the simulation, this time grows 
without limit as the melting temperature is approached, 
consistent with an approach to criticality and a second 
order phase transition. The range of values that was 
explored, however, was not large enough to detect a pos- 
sible power law behavior. These results did not depend 
on the initial conditions used (IC1 or IC2) and on the 
critical temperature used (T c = 0.40725 or T c = 0.4095, 
as there is no exact melting temperature but rather an 
interval [T_,T+]). 
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FIG. 4: Relaxation time as a function of temperature as the 
transition temperature T c is approached from above. There 
are four set of points, corresponding to the two different ini- 
tial conditions (ICl and IC2) indicated in Figure [I] and two 
possible values for T c , determined by the finite interval in 
which T c is found. The time increases without limit as T c is 
approached, consistent with criticality. Error bars are given 

by t min ^max ■ 
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FIG. 2: Upper panel: Pair correlation function at T = 
0.40725 indicates long-range translational order, and at T = 
0.41000, absence of it. Lower panel, orientational correlation 
function at T — 0.40725 indicates long range order, and at 
T = 0.41000, lack of it. 
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Behavior of elastic constants is consistent with 
KTHNY theory 



Another way to determine the nature of the phase tran- 
sition at p — and whether it is consistent with theoreti- 
cal predictions consist on monitoring the behavior of the 
elastic constants of the system. The elastic response of 
an isotropic, homogeneous, continuum solid is character- 
ized by two Lame coefficients A and (i. They appear in 
the compliance tensor as 
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FIG. 3: Final enthalpy as a function of temperature. It 
presents a significant increase at T c — 0.40815 ±0.00090. The 
low temperature phase is a solid, as evidenced by the behav- 
ior of pair and orientational correlation functions. The high 
temperature phase is a liquid (Figure 



tij is the strain and Oki the stress. There are several 
possible ways of extracting A and fi from this tensor. In 
the continuum theory they are of course equivalent. But 
in a numerical calculation involving a finite number of 
atoms, this will no longer necessarily be the case. They 
should coincide, however, within numerical accuracy and 
error bars (see below). The following relations hold for 
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the combination of Lame coefficients: 
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It is a significant prediction of KTHNY theory that K 
approaches a universal value as the critical temperature 
T c is approached from below: 
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where b is the Burgers vector of the dislocations (given 
by the lattice constant at zero temperature). K vanishes 
above T c . 

We have computed the strain of our system following 
Ray and RahmarPlifll, 



e =* [(h^Vhh^-l] 



(8) 



where h is a two-by-two matrix whose column vectors 
define the simulation box, hji is a reference box (here 
taken as the time average of the simulation box) and I 
is the identity matrix. The compliance tensor is given in 
terms of strain fluctuations through 



Sijki = PV R {{eije k i) - (eij)(e k i)) 



(9) 



where Vr is the volume of the reference box. 

The compliances were calculated from simulations of 
2.5 x 10 6 time steps, after 5 x 10 5 equilibration steps. 
The err or ba rs were estimated by performing blocking 
averaged 33 ^ 34 ! and then propagating the error in equation 
The simulation data was divided into 5 data blocks. 
The value of the Burgers vector was estimated as the lat- 
tice constant that minimizes the energy of a triangular 
crystal, bo = 1.11145cr. Figure [5] shows the ratios K\/K c , 
K2 / K c and K3 / K c as a function of temperature as the 
transition is approached from below. The computed val- 
ues are consistent with the KTHNY theory. 



D. Proliferation of dislocations is consistent with 
KTHNY theory 
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FIG. 5: Elastic constant K\/K c (left-hand-side panel), 
K2 /K c (middle panel) and K3 / K c (right-hand -side panel) as 
a function of temperature. Within numerical accuracy they 
coincide, as they should. Near the transition their value is 
consistent with 1, as predicted by KTHNY theory. 



lattice, Pq = 1 and P k = for n / 6. A dislocation 
is characterized by two neighboring sites, one with five, 
and the other with seven, neighbors. The KTNHY the- 
ory predicts that the loss of long range translational order 
is due to the proliferation, and subsequent unbinding, of 
thermally generated dislocation pairs. Such pairs will be 
characterized then by clusters of four sites, two of them 
with five, and two of them with seven, neighbors. Figure 
[6] shows that across the solid-to-liquid transition there is 
a significant decrease in Pq, and a corresponding increase 
in P 5 and P7, consistent with the KTHNY theory. 
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FIG. 6: Pk, the fractional number of sites with k neighbors, 
as a function of temperature for k = 4, 5, 6, 7, 8. There is a sig- 
nificant decrease in P% across the solid-liquid transition with 
a corresponding increase in P5 and P7 , consistent with a tran- 
sition driven by the unbinding of dislocation pairs. 



Finally, and to fully characterize the transition ob- 
served at p = 0, we monitored the number of disloca- 
tions as a function of temperature, counting the number 
of neighbors for each point using a Voronoi construction 
after equilibrium had been reached. Figure [6] shows P k , 
the fractional number of sites with k neighbors, as a func- 
tion of temperature. Of course, for a perfect triangular 



Figure [7] also provides a visual illustration of the num- 
ber of dislocations, monitored with the number of sites 
having 5 or 7 nearest neighbors, within the simulation 
box for different temperatures. Their proliferation is ap- 
parent, consistent with the KTNHY theory. 
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FIG. 7: Dislocation dipole population as a function of tem- 
perature, monitored with the sites having 5 (red) and 7 (blue) 
neighbors. Upper left panel, T = 0.400; upper rigth panel, 
T = 0.405; lower left panel, T = 0.40725, all three temper- 
ature values below T c . Lower rigth panel, T = 0.4095 just 
above T c . The number of dislocation dipoles steadily increases 
as the transition is approached. Sites with 8 neighbors (green) 
are also indicated. 



IV. DISCUSSION: THE ROLE OF EXTERNAL 
PRESSURE 

An hcxatic phase has been observed in simulations pre- 
sented by Chen et alP^ and also reproduced here. How- 
ever, this phase is transient, not in equilibrium, and oc- 
curs when the Lennard- Jones system is subjected to a 
significant external pressure of 20. The use of a constant 
pressure ensemble molecular dynamics ensure the exis- 
tence of a homogeneous phase. However, we have been 
unable to observe the hexatic phase at a constant van- 
ishing external pressure. Is there a reason, within the 
KTHNY theory, not to observe an hexatic phase at zero 
external pressure with a finite number of particles? 

The hexatic phase arises^ after (i.e., at a higher tem- 
perature) a triangular lattice undergoes a dislocation un- 
binding transition but before (i.e., at a lower temper- 
ature) a disclination unbinding transition occurs. The 
latter is possible because the plasma of free dislocations 
screens the disclination-disclination interaction, allowing 
a transition much like the one originally considered by 
Kosterlitz and Thouless". As mentioned in the Intro- 
duction, it is critical, in a numerical simulation, to have 
several length scales available, since the KTHNY mech- 
anism involves the interaction among defects of different 
sizes, at many scales. This is in addition to the fact that 
the theoretical analysis is carried out in the thermody- 
namic limit. We now argue that at least 10 6 particles are 
needed in simulations, in order to have three decades in 
length scales. 

The physics of dislocation unbinding changes consid- 



erably when an external pressure p is included. Indeed, 
the energy U of a dislocation dipole with Burgers vector 
b, whose components are separated by R in this case is^ 

b 2 K " 

log | — ) + C - - cos 2(9 + PbRamO 
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where r is the dislocation core size; C, roughly, deter- 
mines the core energy (i.e., the minimum energy needed 
to generate a dipole); R — \R\\ b = \b\ and 9 is the an- 
gle between R and b. Clearly, when sin# < the last 
term on the right hand side turns the dislocation-dipole 
unbinding process into a thermally activated one, with 
an activation energy Ua (taking 9 = —ir/2, the most 
favorable case, for illustration purposes) 
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Consequently, and given the logarithmic dependence on 
the ratio of external pressure to elastic constant, even a 
modest value of external pressure P, compared to K, will 
give values for the activation energy in the same ball park 
as the chemical potential for the dislocation dipole, and a 
proliferation of isolated dislocations will ensue. Thus, it 
will not be surprising in those circumstances to observe 
an hexatic phase. However, Eqn. (10) shows that, at 



any given non vanishing pressure there will be a finite 
rate of dislocation generation, driving the system away 
from equilibrium. A quantitative study of this interesting 
phenomenon is outside the scope of the present paper. 

In the absence of external pressure, disclination pairs 
above the dislocation unbinding transition interact via 
an energy that depends on the logarithm of their mu- 
tual distance, with a coupling (called Ka by Nelson and 
HalperirP) that is finite due to the screening effect of the 
free dislocations. A second transition towards the liquid 
state thus occurs because of two distinct screenings: free 
dislocations screen the interaction between disclination 
pairs to an effective logarithmic interaction, and then 
this interaction is rcnormalized because of the interac- 
tion between disclination pairs at different length scales. 
So, three length scales should be needed for this scale de- 
pendent interaction among disclination pairs to become 
operative. In addition, a further length scale would seem 
to be necessary in order to have enough free dislocations 
in between a disclination pair for their interaction to be 
effectively screened. According to this reasoning, at least 
~ 10 6 particles would be needed to observe an hexatic 
phase as an equilibrium phase at zero external pressure 
in a numerical simulation. 



V. CONCLUSIONS 

The KTNHY theory of melting in two dimensions'^ 
involves the interaction among dislocation dipoles whose 
sizes span different length scales. Thus, a numerical sim- 
ulation that aims at verifying the theory should involve 
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enough different length scales for this interaction to be 
possible. In two dimensions, 10 particles thus appears as 
an a absolute minimum. Increasing this number should 
improve the statistics. By the same token, 10 6 particles 
should be a minimum number to capture this type of 
effect in three dimensions. Our simulation has been car- 
ried at constant (vanishing) external pressure in order 
to prevent phase coexistence. We find a solid-to-liquid 
transition, and the behavior of the solid phase as the 
transition temperature is approached is consistent with 
the predictions of KTNHY. The behavior of enthalpy as a 
function of temperature is less conclusive: it changes sig- 
nificantly across a narrow temperature range AT, from a 
solid low temperature phase to a liquid high temperature 
phase. The behavior within AT could not be resolved be- 
cause of the limited time-scale that can be reached with 



simulations. There could be an abrupt discontinuity, as 
in a first order transition, or there could be a smooth 
change, including a temperature range with an hexatic 
phase. Above T c , the relaxation time increases as T c is 
approached, consistent with criticality. 
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